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ABSTRACT 

Accurate photometric and kinematic modelling of disc galaxies requires the inclusion 
of radiative transfer models. Due to the complexity of the radiative transfer equation 
(RTE), sophisticated techniques are required. Various techniques have been employed 
for the attenuation in disc galaxies, but a quantitative comparison of them is difficult, 
because of the differing assumptions, approximations and accuracy requirements which 
are adopted in the literature. In this paper, we present an unbiased comparison of four 
methods to solve the RTE, in terms of accuracy, efficiency and flexibility. We apply 
them all on one problem that can serve as a first approximation of large portions of 
disc galaxies: a one-dimensional plane-parallel geometry, with both absorption and 
multiple scattering taken into account, with an arbitrary vertical distributions of stars 
and dust and an arbitrary angular redistribution of the scattering. We find that the 
spherical harmonics method is by far the most efficient way to solve the RTE, whereas 
both Monte Carlo simulations and the iteration method, which are straightforward to 
extend to more complex geometries, have a cost which is about 170 times larger. 
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1 INTRODUCTION 

In order to study the structure of galaxies, it is necessary to 
obtain their intrinsic three-dimensional light distribution, 
i.e. to deproject their two-dimensional image. This depro- 
jection is complicated by the effects of interstellar dust, 
which change the radiation field on its way through the 
galaxy by absorption, emission and scattering. The amount 
of interstellar dust present in disc galaxies has been the 
subject of debate for the last few years. The most widely 
supported idea is that spiral galaxies are moderately opti- 
cally thick in their central regions (a face-on optical depth 
of order unity in the V band), and optically thin in the 
outer regions. This idea is supported by statistical stud- 
ies (Giovanelli et al. 1994, Boselli & Gavazzi 1994, Buat 
& Burgarella 1998), studies of the extinction of galaxies in 
overlapping pairs (White & Keel 1992, Berlind et al. 1997, 
Ronnback & Shaver 1997, Pizagno & Rix 1998, White et al. 
2000), and detailed modelling of the extinction in individual 
galaxies (Jansen et al. 1994, Ohta & Kodaira 1995, Xilouris 
et al. 1997, 1998, 1999, Kuchinski et al. 1998). However, 
there is still no consensus, and other authors claim that spi- 
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ral galaxies are optically thick all over their disc (Valentijn 
1990, 1994, Burnstein et al. 1991, James & Puxley 1993, 
Peletier et al. 1995). Moreover other arguments complicate 
the discussion: it is a gross simplification to talk about the 
opacity of spiral galaxies, because there can be a large dif- 
ference in opacity between arm and interarm regions (White 
et al. 2000). For a detailed overview of the opacity of disc 
galaxies we refer to Davies & Burnstein (1995) and Kuchin- 
ski et al. (1998). 

Many authors have demonstrated that dust attenuation 
(i.e. the combined effect of absorption and scattering) has 
considerable effects on photometric properties such as mag- 
nitudes, colors and scalelengths (Witt et al. 1992, Byun et al. 
1994, Corradi et al. 1996). These effects are of a complicated 
nature, and are not simply proportional to the optical depth 
of the galaxy. Also kinematic studies can be complicated by 
dust attenuation because the projected kinematics will be 
biased towards the motions of stars on the near side of the 
line-of-sight (Bosma et al. 1992, Matthews & Wood 2000, 
Baes & Dejonghe 2001a). Both photometric and kinematic 
studies of disc galaxies hence require sophisticated deprojec- 
tion techniques which take dust attenuation into account. 
The only means to do this properly is by constructing ra- 
diative transfer models. 

The radiative transfer problem is a well-defined prob- 
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lem. It is described by the radiative transfer equation (RTE) , 
which requires as input a precise knowledge of the optical 
properties of the dust grains, the opacity (the total amount 
and the spatial distribution of the dust) and the emissivity 
(the spatial and energy distribution of the stars). If a con- 
venient algorithm is found to solve the equation, the output 
is the projected image on the sky, or the surface bright- 
ness. But usually the problem is the reverse. Given the ob- 
served image on the plane of the sky, can we recover the 
three-dimensional distribution of stars and dust ? The most 
straightforward way to solve this problem is to construct a 
large set of models with various parameters, solve the trans- 
fer equation for each of these models, and then determine 
the parameters such that the obtained solution fits the ob- 
servations. It is clear that efficient algorithms to solve the 
RTE are necessary. 

Unfortunately, the RTE is a fairly complicated equa- 
tion, and it is not straightforward to solve it, unless some 
simplifying assumptions are made. The most widely adopted 
way to simplify the RTE is an approximation or down- 
right neglect of the scattering by dust grains (Guiderdoni 
& Rocca-Volmerange 1987, Disney et al. 1989, Calzetti et 
al. 1994, Ohta & Kodaira 1995). Nevertheless many efforts 
have been made to develop methods to solve the RTE ex- 
actly. Four of the methods that have been explored are 

• the spherical harmonics method, where all terms in the 
RTE are expanded into a series of spherical harmonics, such 
that the RTE is replaced by a set of ordinary differential 
equations. 

• the discretization method, a technique borrowed from 
the stellar atmospheres theory. In this method integrals are 
replaced by sums and the differentials by finite differences, 
resulting in a set of vector equations which can be solved 
iteratively. 

• the iteration method, where intensity is expanded in a 
series of partial intensities. Each of the partial intensities 
obeys its own RTE, which can be solved iteratively. 

• Monte Carlo simulations, where the transfer of photons 
through the gala^cy is investigated by examining the individ- 
ual paths of a large number of photons. 

Each of these methods will have its advantages and disad- 
vantages. On the one hand the different approaches have 
different physical backgrounds, and the selection of a cer- 
tain algorithm can be useful to gain physical insight into the 
problem. For example, the spherical harmonics method gives 
direct access to the moments of the radiation field, whereas 
the iteration method yields direct information about the im- 
portance of multiple scattering. On the other hand, com- 
putational and practical considerations can be motives to 
prefer one method above another. However, these consider- 
ations are often based on general truths without a quan- 
tification. For example, it is generally accepted that Monte 
Carlo methods are costly, but how costly in respect with 
other methods ? 

It is important to be able to estimate how a certain 
solution method scores in terms of accuracy, efficiency and 
ffexibility. According to one's specific interest (e.g. providing 
a tool for statistical studies of attenuation, or constructing 
detailed radiative transfer models for individual galaxies), 
one of these properties may be more important than another. 
Knowing the relative performance of the different methods 



then allows to select the most suitable candidate for the 
problem. 

The aim of this paper is to present an unbiased com- 
parison of the four methods considered. All of them have 
already been applied to construct radiative transfer mod- 
els in order to investigate the attenuation in disc galaxies 
(references will be given). However, the assumptions con- 
cerning the geometry and the dust properties made by the 
various authors often differ significantly. This makes a com- 
parison of the various adopted methods in the literature in 
terms of accuracy, efficiency and flexibility next to impossi- 
ble. Therefore we adapt them such that they can be used to 
solve the same radiative transfer problem. We restrict our- 
selves to a one-dimensional plane-parallel geometry, however 
with absorption and multiple scattering properly accounted 
for. Moreover, all of these models allow an arbitrary vertical 
distribution of stars and dust, and an arbitrary angular re- 
distribution function (ARF), such that they can serve as a 
first approximation to model large portions of galactic discs. 

In Section 2 we discuss the radiative transfer problem 
and present a disc galaxy model to illustrate the solution 
techniques, which are presented in the next four sections. 
In Section 3 the RTE is solved by an expansion in spher- 
ical harmonics, in Section 4 we describe the discretization 
technique, in Section 5 the RTE is solved by an iterative 
method, and in Section 6 a Monte Carlo simulation is used. 
In Section 7 we compare the different techniques in terms of 
accuracy and numerical efficiency, and we discuss them in 
the light of their use to construct detailed models for disc 
galaxies. Section 8 sums up. 



2 THE RADIATIVE TRANSFER PROBLEM 

2.1 The RTE in plane-parallel geometry 

In plane-parallel geometry, there are only two independent 
variables necessary to determine a position and a direction 
in the galaxy: the depth z, and ^, the cosine of the an- 
gle between that direction and the face-on direction (plane- 
parallel geometry implies an azimuthal symmetry around 
the face-on direction ^ = 1). The general form of the trans- 
fer equation in plane-parallel geometry can be written as (in 
this paper we omit all explicit reference to the wavelength 
dependences) 

dl 

+ ^lok{z) £^ I{z, m')*(m, m') V- (1) 

Here I{z,^) is the specific intensity of the radiation, n{z) is 
the dust opacity, rjt{z) is the stellar emissivity, ui the dust 
albedo and 'ii{fi,jj,') is the ARF. The last term in equa- 
tion (|^) represents the scattering term, i.e. the net amount 
of photons which are scattered off dust grains into the direc- 
tion fi from all other directions fi' at z (see Appendix A) . We 
assume that the scattering is coherent, i.e. there is no energy 
redistribution of the scattered light. The thermal emission 
of the absorbed light by the dust grains is neglected, because 
this occurs at far-infrared wavelengths, whereas we will fo- 
cus on shorter wavelengths. Furthermore we assume that the 
dust grains have the same properties over the galaxy. This 
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means that the albedo and the redistribution function are in- 
dependent of position, and that both the emissivity and the 
opacity are separable functions of position and wavelength. 

We replace the height z above the galactic plane by the 
optical depth r, defined as 



f oc 

t{z) — / K{z')dz' , 

J z 



and the RTE then becomes 



(2) 



ti^ir, ^) = J(r, ^) - S,{t) - I J'^ I{r, m') V, (3) 

where 5'*(t) = r7«(r)/K(r) is the stellar source function. This 
equation is to be solved for < r ^ ro, the total optical 
depth of the slab, 



TO = 



K(z)dz 



(4) 



The boundary conditions for our radiative transfer problem 
are 



J(0,-^t) = J(ro,At) = for^i>0. 



(5) 



which means that there is no incident radiation on either 
side of the galaxy. The RTE together with the boundary 
conditions (fel) allow us to calculate the radiation field at any 
position and into any direction in the galaxy. In the first 
place however, we are interested in the radiation field /(0,pi) 
that leaves the galaxy at r = into a certain direction 
H, and more precisely in the fraction of this intensity that 
is attenuated. Because the RTE is a linear equation, this 
fraction will be independent of the total amount of stellar 
emission. We can thus choose the normalization of the stellar 
emissivity (or the source function) arbitrarily. We take 



rj, (2) dz ■ 



S.(r)dr : 



(6) 



which means that the intensity that leaves the galaxy in the 
absence of dust in the face-on direction equals 1. In another 
direction jj,, the dust-free intensity is then simply 1/fj, and 
the attenuation (in magnitudes) is 



A(/x) = -2.51og[p/(0,M)]. 



2.2 A disc galaxy model 



(7) 



In order to solve the problem we still need to characterize 
the parameters and functions that appear in (|l|), i.e. the 
optical properties of the dust and the spatial distribution 
of stars and dust have to be specified. Because the present 
paper focuses on the comparison of different methods to 
solve the RTE rather than on an investigation of the impact 
of dust attenuation on disc galaxies, we restrict ourselves 
here to a single galaxy model to illustrate our results. For 
a thorough study of the attenuation curve as a function of 
various parameters determining the distribution of stars and 
dust, we refer to the second paper in this series (Baes & 
Dejonghe 2001b). 

The vertical distribution of stars in disc galaxies is still 
a matter of debate. Stars were first believed to be isother- 
mally distributed in the z-direction (van der Kruit & Searle 
1981), and later on an exponential behavior was proposed 
(Wainscoat et al. 1989). Nowadays it is believed that the 



Table 1 . The adopted values for the optical properties of the dust 
grains. Tabulated are the optical depth r relative to the V-band 
value, the scattering albedo u> and the asymmetry parameter g. 



band 


A (fira) 


r 


UJ 


9 


U 


360 


1.60 


0.57 


0.49 


B 


440 


1.32 


0.57 


0.48 


V 


550 


1.00 


0.57 


0.44 


R 


700 


0.73 


0.54 


0.37 


I 


850 


0.47 


0.51 


0.31 



true distribution lies in between these two profiles (Pohlen 
et al. 2000, Schwarzkopf & Dettmar 2000). We adopt an 
exponential profile, 
1 



\z\/h 



2h 



(8) 



where h is the scaleheight, which satisfies the normalization 
condition (^). For the dust distribution we also assume an 
exponential distribution. In a normal disc galaxy, the inter- 
stellar matter will sink down to the central plane and form 
an obscuring layer which is narrower than the stellar layer. 
Detailed radiative transfer models of edge-on galaxies in- 
dicate that the scaleheight of the dust is about half of the 
stellar one (Xilouris et al. 1999). As opacity function we then 
obtain 



.(2) = -e- 



(9) 



The normalization of (|9|) is such that 1^ is satisfied. For the 
total optical depth we adopt a moderate value tv ~ 1, which 
seems to be appropriate for disc galaxies, as we argued in 
the introduction. 

For 'i/{jJ., jJ.') we adopt the Henyey-Greenstein ARF (see 
Appendix A), appropriate to describe anisotropic scattering. 
It is a one-parameter function with a (wavelength depen- 
dent) parameter g, which is called the asymmetry parame- 
ter and is the average of the cosine of the scattering angle. 
Values for the optical properties of the dust grains (the wave- 
length dependence of the opacity k, the albedo uj and the 
asymmetry parameter g) are those calculated by Maccioni 
& Perinotto (1994) and displayed in Di Bartolomeo et al. 
(1995). The adopted values of the U, B, V , R and / bands 
are tabulated in Table ^ 

We wish to stress again that this galaxy model is only 
illustrative, and that the methods discussed in this paper are 
applicable to any galaxy model of the kind we consider here, 
i.e. with arbitrary vertical distribution of stars and dust and 
with arbitrary ARF. 



3 THE SPHERICAL HARMONICS METHOD 

One of the most popular techniques to solve the radiative 
transfer equation is a method which uses an expansion in 
spherical harmonics. Many papers have been written which 
describe this method in detail, e.g. Davison (1957), Flannery 
et al. (1980), Roberge (1983), Di Bartolomeo et al. (1995). 
Our approach is based on Roberge (1983), whose method 
can account for any vertical distribution of stars and dust, 
and uses the adaptations of Di Bartolomeo et al. (1995) in 
order to make the matrix in the eigenvalue problem sym- 
metric. 
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In plane-parallel symmetry, spherical harmonics reduce 
to Legendre polynomials, and we expand the intensity and 
the redistribution function in a series of Legendre polyno- 
mials, 



I{r,fi)=J2(^l + l)fiiT)n{fi) 



(10) 



(11) 



galaxy model consisting of a homogeneous mixture of stars 
and dust. Their technique is based on a discretization of the 
RTE on a fixed mesh of points, and difference equations re- 
place the differential equations. We extend the technique of 
Bruzual et al. (1988) such that it can account for any ver- 
tical distribution of stars and dust, i.e. any source function 



S«(r). As we will discuss in section 4.2, we will have to make 



a distinction between source functions that remain finite ev- 
erywhere and source functions that diverge in places. 



The coefficients /i(r) are unknown functions, whereas the 
coefficients ai are known. For example, in the case of 
Henyey-Greenstein scattering, cr; = p', with g the asymme- 
try parameter (Appendix A). Inserting the expansions (|lo|) 
and (|ll|) in (|^, and using the recurrence relations for Leg- 
endre polynomials, we find 

oo 

[lfU{r) + + l)//+i(r) - hif,{T)\Pi{ii) = -S',(r),(12) 

where we set hi = [21 + 1)(1 — uoai). Defining the functions 
ipi and gi as 

Mr) = ^ifi{T) (13) 
g,(T) = -^5i,o, (14) 



we find an infinite set of linear first-order differential equa- 
tions 

' V^;_,(r) + 4±i=^;+i(r)=^,(r) + g,(r). (15) 



\/hi+ihi 



We adopt the so-called Pl approximation (consists in as- 
suming iPi{t) = for I > L, for a certain value of L, L odd) 
in order to turn this infinite set into a finite one, which we 
can write as a vector equation 



AiP'{t) = ilj(r)+g{T). 



(16) 



The RTE is thus reduced to a set of ordinary differential 
equations, for which the system matrix ^ is a non-singular, 
symmetric, tridiagonal matrix with fixed coefficients. Such 
a problem is best solved by diagonalization of the matrix, 
i.e. an eigenvalue problem. The procedure to obtain the ex- 
pressions for ij^iij), are nicely described in Roberge (1983). 

We note in particular that for the Pl solution, the two 
boundary conditions (^ can only be satisfied for (L + l)/2 
directions /ii, which are the positive zeros of the Legendre 
polynomial of order L The intensity 



97-1-1 

/(r,M)=^^^K^)^'KM). 



/hi 



(17) 



at other directions can be obtained by a cubic spline inter- 
polation. In general we found that the Pl solution converges 
very fast, and that the boundary conditions were met when 
L ^ 15, which we adopted in our calculations. This value is 
slightly larger than that used by Di Bartolomeo et al. (1995). 



4 THE DISCRETIZATION METHOD 

One of the first major efforts to investigate the effects of 
absorption and multiple scattering in disc galaxies, is the 
paper by Bruzual et al. (1988). They solve the RTE for a 



4.1 Finite source functions 

If the source function S,{t) remains finite for all values of 
r, we can easily generalize the procedure of Bruzual et al. 
(1988). Starting from the transfer equation ^ we introduce 
the even and odd fields 

^(r,/.) = i[/(T,M)+/(T,~M)] (18) 
v{r,^) = ^[I(r,f^)-IiT,^^i)]. (19) 

The transfer equation can then be written as 



^('^' ^) = ^('^' A*) - ^ / ^{t, a*')*°(m, ^J■')<i^^' 
OT Jo 



(20) 



dv , 



A* ^("^'Z^) = w(r,/i) - 5'*(r) j u{t,ij, )*''(/i,Ai )d/i , 



where 

*=(M,M') = f [*(m,m') + *(m,V)] 

*°(m,m') = H*(m,m')-*(m,V)] 

The new boundary conditions are 



(21) 

(22) 
(23) 

(24) 
(25) 



These equations are analogous to the ones Bruzual et al. 
(1988) use for their homogeneous slab, with the exception 
that the constant source function for the homogeneous slab 
is replaced by a non-constant, but finite source function 
S*{t). a similar discretization procedure can be followed. 
We introduce a grid of A'^ -I- 1 mesh points r^, uniformly 
spaced in optical depth space. In between these points we 
introduce a second grid, denoted by half-integer numbers 
Ti+i/2- A derivative evaluated in a half-integer grid point 
Ti+i/2 becomes a difference evaluated in the nearest inte- 
ger mesh points n and r^+i, and analogous for the integer 
mesh points. At the boundaries tq and rjv the derivatives are 
expressed by means of the boundary conditions. The inte- 
grals over fj, are approximated by M-point Gauss-Legendre 
quadrature, and we hence introduce a mesh of M points ^j, 
being the roots of the Mth order Legendre polynomial. The 
RTEs are then replaced by a set of linear vector equations, 
where the {2N + 1)M unknowns are the even and odd fields 
Uij and Vi^i/2j evaluated at the integer and half-integer 
mesh points respectively. These equations are solved recur- 
sively using (a simplification of) the elimination scheme of 
Milkey et al. (1975). 

Optimal values for M and N are hard to determine. 
On the one hand these values should not be too high, be- 
cause of computational limitations. Indeed, the elimination 
scheme consists of two loops, and every step in the first loop 
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consists of 2 matrix inversions, where the order of the ma- 
trices is M, the number of angle quadrature points. This 
means that totally 2N + 1 matrices of order M need to be 
inverted, which is a numerically costly process. Moreover, 
the memory requirement is high, because the 2N + 1 matri- 
ces calculated in the first loop need to be stored, for use in 
the second loop. On the other hand, the values for M and 
should not be taken too low, such that the approximation 
of differentials by differences and integrals by quadrature 
sums is acceptable. Typically M — 10, while typical values 
for TV are dependent on the total optical depth tq and should 
be chosen such that At never exceeds 0.01. For our galaxy 
model this meant > 100. 



4.2 Diverging source functions 

The above procedure cannot be used anymore when the 
source function diverges at the boundaries t — and/or 
T — To, because the source function needs to be evaluated 
in these end points. Such source functions are realistic, be- 
cause they correspond to all distributions where the dust 
distribution is narrower than the stellar distribution, which 
is observed in most galaxies. In particular, the galaxy model 
described in Section 2.2 has a diverging source function. 

One obvious solution would be to apply a cutoff, but the 
method then becomes very dependent on the actual position 
of the cutoff. Moreover such a solution is inaccurate, because 
it involves a lot of matrix inversions, and these are difficult 
operations if the matrices contain elements of very difi'erent 
magnitudes. 

If the source function diverges, we don't use the opti- 
cal depth coordinates, but we define a similar independent 
variable ^, 



77,(2') dz'. 



(26) 



which assumes values between and 1, because of the nor- 
malization of the emissivity The RTE (jl|) becomes 



^LoR,iO |'^/(C,M')*(M,M')dM', (27) 



where is the reciprocal of the stellar source function. 



R*iO 



1 _ <0 



(28) 



Starting from this form of the RTE, we can repeat the 
discretization procedure from the previous paragraph, i.e. 
we construct a uniform mesh and an intermediate mesh 



e 



i+1/2 



for the finite differences, and a mesh 



for the 



Legendre-Gauss quadrature, and we obtain a set of vec- 



tor equations for the unknowns Uij and v. 



which are 



solved by the elimination scheme of Milkey et al. (1975). 

The arguments leading to the choice of A'' and M are 
the same as in the previous paragraph. 



5 ITERATION METHODS 



5.1 The intensity as a series 

In a very early paper on scattering, Henyey (1937) shows 
that the RTE (^ can be solved iteratively by writing the in- 
tensity as a summation of partial intensities 7„, which repre- 
sent the radiation field consisting of photons that have been 
scattered exactly n times. This technique has been adopted 
by van de Hulst & de Jong (1969) and Xu & Helou (1996) 
to solve the RTE in plane-parallel geometry. Specifically in 
our case, we write the intensity as 



n=0 

The nth partial intensity satisfies the RTE 
dl„ 

= In{T,y) - Sn{T,y), 

where 

So{t,ii) = S^t) 

S'„(r,/i) = ^ j I„-i{T,n')'^{fi,n')dfi', 

and is subject to the boundary conditions 
/„(0, -n) = In{To,fi) = for ^ > 0. 



The solution of 



5^) can be directly written 
5'„(r',^) exp (- dr' 

r 

S„{t',ii) exp 



dr'. 



(29) 

(30) 

(31) 
(32) 

(33) 

(34) 
(35) 



We calculated the partial intensities and source functions 
on a mesh of size L in optical depth and size M in angle. 
Because the integration over the angle in (^^ always has the 
same integration domain, we chose the M angle points as the 
abscissae for an A/-point Gauss-Legendre quadrature. For 
the calculation of the I„, the integration along the path has 
variable boundaries, and hence a simple quadrature would 
yield bad results near the edges of the mesh. Instead, we used 
a uniform mesh of optical depth points, and for each fixed 
angle ^j, we constructed a cubic spline approximation to 
S„{t, fij), and performed the integration using this function. 
Typical values for L and M are 40 and 10 respectively. 

The number of terms necessary in the summation ( |29[ ) 
in order to obtain an accurate result depends on the wave- 
length, because the nth term in the expansion is propor- 
tional to the nth power of the albedo . Xu & Helou (1996) 
used 20 terms for their sandwich model; we find that in gen- 
eral, about 10 terms are sufficient. For our one-dimensional 
plane-parallel geometry this number of integrations is still 
manageable. However, if the method is extended to more 
complicated geometries, more dimensions will be added, and 
calculation of such a high number of terms becomes very 
time-consuming, such that adaptations of the method are 
desirable. 



5.2 The intensity as a geometric series 

A simplification of the iteration method was introduced by 
Kylafis & Bahcall (1987). They use the approximation that 
the amount of photons that are scattered exactly n-l-1 times 
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Figure 1. The ratio Pn{^i) = l,i(0, a)/Jp i (0. u) at tlie V band 
of two consecutive terms in the series (G9|), as a function of n. 
It is shown for four different inclination angles, for the model 
described in Section 2.2. 



to the amount of photons scattered exactly n times, is con- 
stant, i.e. 



/„+i(r,^) _ In{r,fJ.) 



7„(r, ^) /„_i(r,/i) 



for all n. 



(36) 



If this approximation holds, only two terms in the sum ( |29[ 
need to be calculated, because it can be written as a geo- 
metric series 



3=0 



h (r, fi) 



(37) 



To test the accuracy of the approximation (|3^), we calcu- 
lated several terms in the sum (E9|) for our galaxy model. 
The results are shown in Figure H7where we plot the ratio 
Pn(M) = /n(0, /x)//n-i(0, /^) as a function of n, for four dif- 
ferent values of ^. It is clear that Pn{j-i) is indeed constant for 
n > 3, and that this constant is independent of the angle fi. 
For the first values of n however, the condition (^) is clearly 
not satisfied. Therefore we adapt the strategy of Kylafis & 
Bahcall (1987), and we assume that the condition (^) is 
valid for n > N. The intensity can then be written as 



H'T,!^) = h{T,^j) + Ii{r,fi) + ■ 



+ /jv-i(r, fi)Y^ 

3=0 



lN{r, fi) 



i(^,m) 



(38) 



The parameter A'' thus determines the last of the partial 
intensities that needs to be calculated explicitly. We find 
that already for N = 2 the results are correct with less than 
a one-percent error. The number of integrations that then 
has to be performed for the solution of the RTE is 3LM. 



6 MONTE CARLO SIMULATION 

The last method we considered to solve the RTE is a Monte 
Carlo simulation, which is probably the most widely adopted 
method for radiative transfer problems. The principles of 
this method are described in detail by Cashwell & Ev- 
erett (1959), Mattila (1970), Witt (1977), Yusef-Zadeh et 
al. (1984) and Bianchi et al. (1996). 



6.1 The principles 

The Monte Carlo method basically follows the individual 
path of a very large number A*" of photons through the 
galaxy. At each moment, the fate of a photon on its path 
is determined by a number of quantities such as the free 
path length between two interactions, the nature of the in- 
teraction (scattering or absorption) and the direction change 
during a scattering event. Each of these quantities can be de- 
scribed by a random variable, taken from a particular prob- 
ability density p{x)dx. 

More specifically, in plane-parallel geometry a photon is 
characterized by two variables: the position (or equivalently 
the optical depth r) and the direction fi. To start, the initial 
position Ti and direction /ii are generated randomly as 







(39) 
(40) 



f^^ dr 

X2^J — or ^ii = 2X2-1, 

where Xi and X2 are uniform deviates. Next we generate a 
random free path length £ (also in optical depth units) by 
setting 



X3 



e ^da; or £ -- 



HI-X3) 



(41) 



This randomly determined path length is to be compared 
with the maximal free path length L of the photon in con- 
sideration. 



L = 



Tl 
Ml 



t q — ri 
Ml 



if /ii > 
if /ii < 0. 



(42) 



If £ > L, the photon will leave the galaxy, and its direction 
Hi is recorded. If £ < L, the photon will interact with a dust 
grain. The nature of this interaction can be determined by 
chosing a uniform deviate X4 and setting 



interaction 



scattering 
absorption 



ii X4, < uj 
if X4 > 00- 



(43) 



If the interaction is an absorption, the photon disappears 
and will not contribute to the final intensity. If the inter- 
action is a scattering, the photon will have a new position 
and a new direction. Given the free path length £ and the 
original position t\ and direction ^1 , the new position of the 
photon is 



7-2 = Tl — H\£, 



(44) 



whereas the new direction fj,2 is determined by the uniform 
deviate 



X. 



>l/(/ii, a;) da;. 



(45) 



This procedure can now be repeated until the photon is ei- 
ther absorbed or leaves the galaxy. In the latter case it will 
contribute to the observed intensity and its direction must 
be recorded. After a sufficiently large number of such ex- 
periments, an M-binned histogram of the emerging angular 
distribution can be constructed. Due to the planar symme- 
try of our galaxy models, photons leaving the galaxy at the 
back side can be added to those leaving at the front side. 
The ratio of the number of photons to the number of photons 
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that would be in the bin if there were no dust attenuation 
(i.e. the distribution of the initial directions), then yields the 
attenuation in that bin. The final attenuation curve A{fi) is 
determined by fitting a smooth curve to these results. 



6.2 Optimizing the routine 

There are two simple ways in which the above described rou- 
tine can be optimized, i.e. adapted such that better statistics 
can be obtained with less photons. 

The first adaptation is to avoid that photons disappear 
from the radiation field by absorption. This can be done by 
assigning a weight to each photon. At each interaction the 
probability that the photon will be scattered is uj, whereas 
the probability for absorption is 1 — a;. Instead of simulat- 
ing absorption as described above, we let the photon scatter 
and we reduce its weight by a factor lu. This way, no pho- 
tons disappear from the radiation field, and each photon will 
eventually leave the galaxy. At that point, both direction 
and weight are recorded. The final histogram is obtained by 
counting the number of photons in each bin, weighted by 
their individual weight. 

A second adaptation is the concept of forced first scat- 
tering, by which each photon is forced to be scattered at 
least once before it leaves the galaxy. Given the initial max- 
imal free path length L(r, /i) of the photon, the probability 
that the photon leaves the galaxy is exp(— L). When the to- 
tal optical depth of the galaxy is low, many photons leave 
the galaxy without interaction, such that a large number 
of photons is necessary to obtain reliable statistics of the 
scattered radiation. Hence, instead of allowing that a pho- 
ton can escape from the galaxy, we split it. A fraction with 
weight w = exp(— L) leaves the galaxy and will be classified. 
The other fraction with weight w = — exp(— L)] is forced 
to scatter before it leaves the galaxy. Therefore a free path 
length £ has to be generated such that £ < L. This can be 
done by replacing (klj) with 

X-i = J e~^d2: ^ J e~^da; 

or £ = - In [l - X;i (l - e"^)] . (46) 

After this forced first scattering, the following scatterings 
are as described above. 

These two concepts reduce the number of photons, nec- 
essary to obtain reliable statistics, significantly. Typical val- 
ues of A'^ and M we use in our calculations are N = 10^ and 
M = 100, such that the mean number of photons in each 
bin is around 1000. 



7 DISCUSSION OF THE ADOPTED 
TECHNIQUES 

7.1 Comparison of the numerical results 

Because we have at our disposal four different methods to 
solve the RTE, it is straightforward to check the correct- 
ness of the individual methods by comparing their results. 
In Figure^ we plot the V band attenuation for the model de- 
scribed in Section 2.2, for the four methods considered. We 



Table 2. A check on the accuracy of the four adopted methods 
of solving the RTE. We computed the attenuation A{fi) of the 
model described in Section 2.2, for four angles and for the five 
bands U, B, V, R and I. The four results correspond to the 
spherical harmonics method (sh), the discretization method (di), 
the iteration method (it) and the Monte Carlo simulation (mc). 







U.zo 


U.oU 


0.75 


l.UU 


u 


sh 


1.113 


0.617 


0, 


,356 


0.209 




di 


1.113 


0.617 


0. 


.355 


0.209 




it 


1.119 


0.617 


0. 


.354 


0.208 




mc 


1.108 


0.616 


0. 


.360 


0.213 


B 


sh 


1.010 


0.520 


0, 


,278 


0.149 




di 


1.010 


0.520 


0, 


,277 


0.148 




it 


1.014 


0.519 


0, 


,276 


0.148 




mc 


1.003 


0.524 


0, 


,277 


0.129 


V 


sh 


0.865 


0.398 


0, 


,187 


0.081 




di 


0.866 


0.398 


0, 


,187 


0.081 




it 


0.867 


0.397 


0, 


,186 


0.081 




mc 


0.851 


0.402 


0, 


,188 


0.078 


R 


sh 


0.721 


0.299 


0, 


,123 


0.038 




di 


0.721 


0.298 


0, 


,123 


0.037 




it 


0.722 


0.298 


0, 


,123 


0.038 




mc 


0.710 


0.300 


0, 


,127 


0.029 


I 


sh 


0.537 


0.197 


0, 


,068 


0.007 




di 


0.538 


0.196 


0, 


,067 


0.006 




it 


0.538 


0.196 


0, 


,067 


0.007 




mc 


0.533 


0.195 


0, 


,070 


0.007 



can obviously conclude that the four modelling procedures 
are in complete agreement with each other. 

However, as already mentioned, the modelling tech- 
niques do not reveal the solution at any direction fi. The 
spherical harmonics method only yields A for the (L -I- 1)/2 
positive zeros of the Legendre polynomial of order L + 1. The 
discretization and iteration methods yield A at the abscis- 
sae of the adopted quadrature formulae. The Monte Carlo 
method actually yields a histogram of the attenuation in 
each of the bins into which the interval of possible /x-values 
is divided. In order to find the attenuation in a randomly 
chosen direction fi, we use a spline interpolant for the first 
three methods, and a fitting polynomial for the Monte Carlo 
method. In Table ^ we tabulate the attenuation curve A{^), 
calculated by each of the four methods, for four inclination 
angles fi. This table shows that the differences between the 
different solutions is always of the order AA « 0.01 mag. 
Figure ^ and Table |^ prove that the four methods are accu- 
rate and consistent. 

In principle, we could also check the correctness of 
the methods by comparing our results with those obtained 
by other teams, who adopted similar techniques. Such a 
comparative analysis is conducted by Di Bartolomeo et al. 
(1995), who compared their results for a homogeneous slab 
with those of Guiderdoni & Rocca-Volmerange (1987), Ky- 
lafis & Bahcall (1987), Bruzual et al. (1988), Witt et al. 
(1992), Byun et al. (1994) and Calzetti et al. (1994). The 
discrepancies between the extinction curves are significant 
(e.g. AAu up to 0.3 mag for tv as small as 0.5), and it is 
important to investigate why this is so. Are they due to the 
adopted solution technique or to other causes ? Generally, 
the discrepancies can be the result of differences in 
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Figure 2. Comparison of the V band attenuation curve Av{n) obtained by the four different methods. The results arc only shown for 
the particular values for which the different methods yield the solution directly (see text). 



(i) the solution method. In five of the mentioned papers 
the RTE is solved exactly (i.e. taking absorption and scat- 
tering fully into account) using the discretization, iteration 
or Monte Carlo techniques, whereas Guiderdoni & Rocca- 
Volmerange (1987) and Calzetti et al. (1994) use approxi- 
mate analytical solutions. 

(ii) the grain properties. The fact that different authors 
use different sets of optical properties can introduce substan- 
tial discrepancies in the attenuation curve. These differences 
between the various values for the grain properties can be 
very substantial, as clearly shown in Figures 1 and 2 of Di 
Bartolomeo et al. (1995). 

(iii) the geometry. As Di Bartolomeo et al. (1995) indi- 
cated, it is sometimes impossible to compare results because 
the geometrical configuration used by the various authors 
is not always the same. Witt et al. (1992) use a spherical 
symmetry, Kylafis & BahcaU (1987) and Byun et al. (1994) 
adopt a axisymmetric galaxy model, whereas the other au- 
thors use a plane-parallel homogeneous slab. 

In fact in only two of the seven papers, the RTE is solved 
exactly for a plane-parallel homogeneous slab: Bruzual et 
al. (1988) by using the discretization technique and Di Bar- 
tolomeo et al. (1995) by adopting the expansion in spherical 
harmonics. However, the authors adopt a different set of op- 
tical properties of the dust grains, and the differences AA{ii) 
between the attenuation curves can still be due to the first 
two points mentioned above. 

To test to which degree the adopted technique con- 
tributes to the differences in the extinction curve, we con- 
sider a plane-parallel slab, and solve the RTE twice for each 
of the four methods at our disposal: once with the dust 



grain parameters of Di Bartolomeo et al. (1995) -the ones 
we adopt throughout this paper-, and once with those of 
Bruzual et al. (1988). For a galaxy with total optical depth 
TV = 1, the results are shown in Figure ^ for the U and I 
bands. We find a very good agreement between the results 
obtained by Bruzual et al. (1988) and Di Bartolomeo et 
al. (1995) respectively, and our results with the correspond- 
ing dust parameters. This demonstrates that the attenua- 
tion differences AA are completely due to the differing dust 
properties, and that the spherical harmonics and iteration 
methods are reliable. 



7.2 Computational efficiency 

Besides being accurate, another desirable quality of a solu- 
tion method is the computational efficiency, as we explained 
in the introduction. 

For the spherical harmonics method the only numeri- 
cally costly operations are the calculation of the eigenvalues 
of a matrix of order L + 1, the inversion of such a matrix, 
and L + 1 simple integrations. Given L ~ 15 this cost is 
relatively low. For the discretization method 2A'' -I- 1 matri- 
ces of the order M need to be inverted, with typical values 
Ad = 12 and A'^ = 100. Because the computation time for 
the inversion of a matrix of order Al is proportional to M^ * 
(Press et al. 1989), the numerical cost of the discretization 
method will be considerably higher. The iteration method 
requires 3A4 cubic spline fits, and 3LAI integrations along 
the line-of-sight, with typical values L = 40 and M — 10. 
Although the integrands are well-behaved (the product of 
a smooth source function with an exponential), such that 
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Figure 3. Comparison of our work with that of Bruzual ct al. 
(1988) and Di Bartolomeo et al. (1995). Shown is the attenua- 
tion curve for the U and / bands, for a homogeneous slab with 
Ty = I. The stars and asterixes represent the results obtained by 
Bruzual et al. (1988) and Di Bartolomeo et al. (1995) respectively, 
the curves are our results where we used the corresponding dust 
parameters. 



Table 3. A comparison of the computation time necessary for the 
calculation of the attenuation curve A(^) for a single wavelength. 
The number gives the actual computation time in seconds, the 
second number is the computation time relative to the spherical 
harmonics method. 



spherical harmonics 


0.091 s 


1 


discretization 


2.55 s 


28 


iteration 


15.68 s 


172 


Monte Carlo 


15.25 s 


168 



each individual integration is relatively easy to perform, the 
high number of integrations makes the iteration method nu- 
merically costly. Finally, for the Monte Carlo simulation, the 
only costly operation for each photon trajectory is the gen- 
eration of 2n -)- 3 pseudo-random numbers, where n is the 
number of scatterings for the photon. Because the number 
of photons must be fairly high in order to achieve reliable 
statistics (we use A'^ = 10^), the numerical cost of the Monte 
Carlo method is considerable. 

In Table ^ we tabulate the mean computation time nec- 
essary for the calculation of the attenuation curve for a sin- 
gle wavelength. This little table shows that the iteration 
method and the Monte Carlo simulation have the same ef- 
ficiency, and that the spherical harmonics method is signif- 
icantly more efficient that the other methods. Although for 
this simple one-dimensional plane-parallel geometry the ef- 
ficiency is less important (the computation times are very 
feasible for all methods), it will be important if we want to 
extend these solution to more complex geometries. 



7.3 Extension to more complex geometries 

The implementations described in this paper can handle the 
RTE in a plane-parallel geometry, and although they can 
accomodate any vertical star-dust geometry, they cannot 
model real disc galaxies. More realistic models require a Ught 
distribution that also decreases exponentially in the radial 
direction (Freeman 1970, Saio & Yoshii 1990, Firmani et al. 
1996). Instead of a one- dimensional plane-parallel geometry, 
the RTE then has to be solved in a two-dimensional axisym- 
metric geometry, with a set of four independent coordinates 
instead of a pair. This extra dimension complicates the RTE 
substantially. This observation forces us to think about how 
the techniques described in this paper can be generalized to 
solve the RTE in axisymmetric disc galaxy models. 

The most obvious candidate for extension to axisym- 
metric galaxy models (or any other geometry) is the Monte 
Carlo simulation. This technique has already been applied 
several times to realistic disc galaxy models (e.g. Bianchi 
et al. 1996, Ferrara et al. 1996, de Jong 1996, Wood & 
Jones 1997, Matthews & Wood 2000). Besides being able to 
treat any geometry, it is also sufficiently flexible to treat, for 
example, clumpy dust distributions (Witt & Gordon 1996, 
2000, Bianchi et al. 2000). This flexibility makes Monte 
Carlo simulations probably the most powerful technique for 
solving complicated radiative transfer problems. However, 
their great numerical cost is a disadvantage, and this cost 
will grow strongly if the method is extended to axisymmetric 
geometries. Indeed, given an initial position, direction and 
a pathlength, the next position on the path is directly cal- 
culated in our plane-parallel geometry, but in axisymmetric 
geometries this becomes a very time-consuming operation 
(Bianchi et al. 1996). Therefore it is worth while to investi- 
gate how the other techniques can be generalized. 

Also the iteration technique is easily extendible to more 
complex geometries. And in contrast to the Monte Carlo 
simulation, it can be expected that the numerical cost will 
grow in proportion to the number of dimensions. If meshes 
of order J and K are constructed inorder to account for the 
extra radial and azimuthal dimensions in axisymmetric ge- 
ometry, the entire routine can be expected to take about JK 
times as much computation times. The iteration method has 
been adopted to solve the RTE for axisymmetric disc galax- 
ies in its original form (Vansevicius et al. 1997), but it is 
the modiflcation by Kylafis & Bahcall (1987) that has in- 
creased the efficiency considerably, and turned the method 
into a very practical instrument to investigate dust attenu- 
ation in disc galaxies (Bosma et al. 1992, Byun et al. 1994, 
Misiriotis et al. 2000). In particular, the method has been 
adopted to construct detailed radiative transfer models for 
a set of highly inclined disc galaxies (Xilouris et al. 1997, 
1998, 1999), which have been compared with dust emission 
models (Alton et al. 2000). It is interesting that the results 
of Xilouris et al. (1999) seem at a flrst glance not in corre- 
spondence with the Monte Carlo results of Kuchinski et al. 
(1998), whose computed optical depths are approximately a 
factor 2-3 higher. 

In the previous paragraph, it was shown that, while 
the efficiency of the Monte Carlo and iteration techniques 
is comparable (at least for the one-dimensional geometry), 
the spherical harmonics method is superb in efficiency. It 
is in fact possible to extend the expansion in spherical har- 
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monies to geometries other than plane-paraUel or spherical 
(Davison 1957). However, another way to extend the spher- 
ical harmonics technique to axisymmetric disc galaxies has 
been proposed by Corradi et al. (1996), who assume a local 
plane-parallel geometry along each linc-of-siglit. The RTE 
can then be solved for each line of sight separately, such 
that the computation time will be very small. The limita- 
tion of this method however is that the assumption of a 
local plane-parallel geometry must remain acceptable. For 
(nearly) face-on galaxies this may be the case, because the 
scale at which the galaxy's structure (i.e. the source func- 
tion) changes on the plane of the sky is small compared 
to the mean free path of the photons. For highly inclined 
galaxies however, the mean free path of the photons is large 
compared to the scale of variation of the source function 
on the plane of the sky, such that the assumption of local 
plane-parallel geometry will not be applicable. 

The last described method, the discretization method 
is a typical one-dimensional technique and is least suitable 
for extension to more complex geometries. 



8 CONCLUSION 

Constructing radiative transfer models for disc galaxies re- 
quires a solution of the RTE. This equation is, even for the 
most simple geometries, sufficiently complex, such that so- 
phisticated methods are necessary to solve them. The RTE 
typically takes as input the distribution of stars and dust, 
and the output is the projected image on the sky. How- 
ever, usually the problem is the reverse: the distribution of 
dust and stars has to be derived from the observed inten- 
sity. The most simple way to solve this problem is to find a 
best fitting solution of the RTE in a large parameter space. 
This involves that the RTE has to be solved repeatedly, such 
that efficiency is a very important property of a good solu- 
tion method. For the same reason accuracy is important: 
the eflfects of dust on the attenuation on certain physical 
parameters can be very wealc, such that a high accuracy is 
required to determine the stellar and dust distribution. And 
last but not least, there are many ways to solve the RTE 
for simple one-dimensional geometries, but the modelling of 
real galaxies demands methods which are applicable to ge- 
ometries other than plane-parallel or spherical, in the first 
place axisymmetric. Therefore flexibtUty is necessary. 

Several techniques are adopted in the literature to solve 
the RTE, even for realistic axisymmetric disc galaxy models, 
but the accuracy, efficiency and flexibility of these methods 
has never been properly compared. Nevertheless, it is clear 
that knowledge of these properties is necessary in order to 
select the most suitable method for a particular radiative 
transfer problem. 

In this paper we investigate four different methods to 
solve the RTE in a simple plane-parallel geometry: the spher- 
ical harmonics method, the discretization method, an iter- 
ative method and a Monte Carlo simulation. All of them 
are adapted such that they solve the RTE exactly (i.e. with 
the physical processes of absorption and multiple scattering 
properly taken into account), and that they allow an arbi- 
trary vertical distribution of stars and dust and an arbitrary 
ARE. This way our methods can contribute to understand- 
ing of the transfer of radiation in realistic galactic discs. 



For a galaxy model with realistic vertical structure and 
dust parameters, all four methods yield exactly the same 
results, with differences between the attenuation curves at 
most a few hundredths of a magnitude. They can thus be 
considered as accurate. We also compare our results with 
those obtained by others, and find that the results are in 
agreement with each other. Concerning efficiency, the iter- 
ation method and Monte Carlo method are close, whereas 
the discretization method is 6 times as efficient, and the 
spherical harmonics method even 170 times. Whereas the 
Monte Carlo method can easily be generalized to an arbi- 
trary geometry, we anticipate that the iteration method will 
probably be the most efficient routine in axisymmetric ge- 
ometry. The adaptation of the routine to axial symmetry is 
more or less straightforward, and the efficiency will not suf- 
fer as much from the extra dimensions as the Monte Carlo 
simulation. This issue will be thoroughly investigated and 
presented in a forthcoming paper. 
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APPENDIX A: THE HENYEY-GREENSTEIN 
ANGULAR REDISTRIBUTION FUNCTION 

In general, the phase function "l>(n, n') defines the probabil- 
ity that a photon, that is scattered from a direction n' , will 
obtain a new direction n. If we normalize the phase function 
as 



<I>(n, n' 



dn' 

4tt 



for all n. 



(Al) 




Figure Al. The V band Henyey-Greenstein ARE '^{fj,,fj,') as a 
function of the new direction fi. It is shown for various values of 
the initial direction fi' , ranging between and 1, with intermedi- 
ate step 0.25. For negative values of fi' the ARF is obtained by 
the symmetry relation ^{fi, —fJ.') = ^(— /i, /i'). 



then the amount of photons added to the radiation field 
I{r,n) at a position r into a direction n due to scattering 
is 



ujK.{r) j j I{r,n')^{n,n')— — . 



(A2) 



It can be assumed that the scattering phase function does 
not depend independently on the four variables {n,n') = 
(/i, /i', </)'). We assume that it depends on them only 
through the angle O between the incident and the scattered 
radiation. In our plane-parallel geometry, we have azimuthal 
symmetry, and the scattering term can then be written as 



(A3) 



where the ARF is defined as (van de Hulst & de 

Jong 1969, Bruzual et al. 1988) 



2tv 



$(cose)d<ji'. 



(A4) 



A widely used phase function is the one named after Henyey 
and Greenstein (1941), 



<l>(cos 9) 



{l+g^ - 25cose) 



3/2 ' 



(A5) 



which is an accurate one-parameter family to describe the 
average scattering in galactic dust. Its parameter g is the 
asymmetry parameter and is the mean cosine of the scatter- 
ing angle, g = (cosO). A very useful characteristic of this 
function is that it has a very simple expansion in Legen- 
dre polynomials. Using the generating function for Legendre 
polynomials it is easily shown that 



$(cos e) = ^ {21 + l)g'Pi{cos e). 



(A6) 



such that we immediately have an expression for the angular 
redistribution function (Chandrasekhar 1960, Roberge 1983) 



<fi^^,^,')^Y.('^l + l)g'Pl{^,)Pl{^,') 



(A7) 
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In order to derive a closed expression for '^{fi, jJ.') we use 
cos e = ^i^' + V(l - - A*'^) cos(<;!> - (j)'). (A8) 



and combine this expression with (A4) and (A5), to obtain 

(A9) 



with 
a{fj,,ii') 



1+5 



dt 

(a ± 6 cost) ^''■^ 



= 2|g|V(l-M')(l-M'') 



(AlO) 
(All) 



For all values of the asymmetry parameter, these functions 
obey the relation a > b ^ such that we find (Gradshteyn 
& Ryzhik, 1965, Section 2.575) 



2(1 



n{a — b)\/a + b 



zE 



2b 

a + b 



(A12) 



where E{k) is the complete elliptic integral of the second 
kind. Figure Al shows the Henyey-Greenstein ARF at the 
V band for a few values of the initial direction ^' . 

This paper has been produced using the Royal Astronomical 
Society/Blackwell Science I^T^ style file. 
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